##
## CPS Forthcoming
## "Strategic Logic of Elite Purges in Dictatorships"
##
## Jun Koga Sudduth 
## jun.koga@strath.ac.uk
## November 2016
##
## CPS_Elitepurge_JunSudduth.R
##

library(foreign)
library(MASS)

# update.packages()

setwd("C:/Users/Jun/Dropbox/Elite Purge Project")

dat <- read.table(file="CPS Final/CPS_Elitepurge_Sudduth.txt")
 
dat <- dat[order(dat$ccode, dat$year, dat$leadid2),]

# Recode "startdate" and "enddate"
dat$startdate <- dat$syear*10000 + dat$smonth*100 + dat$sdate
dat$enddate <-  dat$eyear*10000 + dat$emonth*100 + dat$edate_x




##
## Figure 1 
##

b.sim <- read.table("Data/VCVmatrix2_01012015.csv", header=TRUE, sep=",")  # this is a 10000 draws of simulated coefficients.
dim(b.sim) # 10000 * 10
 
f <- coupexit~ Bolivia +  IS_WAR +  chrgdpch  + lrgdpch + logmilex  + mildic + royal + party + failedcoup2_1  + coupentry2*log(tenure+1)
model1 <- glm(f, family=binomial(link = "logit"), data=dat)


#b.sim <- mvrnorm(10000,coefficients(model1),vcov(model1)) # 1000*19 matrix of simulated cofficients.
#dim(b.sim) # [1] 1000    9

coefficients(model1)
tenure.range <- seq(0, 12, 0.25)

# Coup-entry leader

Z <- cbind(1,  median(dat$Bolivia, na.rm=TRUE), 0,  mean(dat$chrgdpch, na.rm=TRUE),  mean(dat$lrgdpch, na.rm=TRUE), mean(dat$logmilex, na.rm=TRUE),  median(dat$mildic, na.rm=TRUE),
median(dat$royal, na.rm=TRUE), median(dat$party, na.rm=TRUE),0,  1,  log(tenure.range +1), 1*log(tenure.range +1))
coefficients(model1)
Zb <- Z%*%coefficients(model1)
P <- exp(Zb)/(1+exp(Zb))      # 201*1 vector of P

Zb.sim <- Z %*% t(b.sim)
P.sim <- exp(Zb.sim)/(1+exp(Zb.sim))   # 201*1000 simulated probabilities.
P.sim.graph <- t(apply(P.sim,1,sort)) # 201*1000 simulated probabilities, sorted to identify the confidence bounds.


# No coup-entry leader

W <- cbind(1,    median(dat$Bolivia, na.rm=TRUE), 0,  mean(dat$chrgdpch, na.rm=TRUE),   mean(dat$lrgdpch, na.rm=TRUE), mean(dat$logmilex, na.rm=TRUE),  median(dat$mildic, na.rm=TRUE),
median(dat$royal, na.rm=TRUE), median(dat$party, na.rm=TRUE),0,  0,  log(tenure.range +1), 0*log(tenure.range +1))

Wb <- W %*%coefficients(model1)

P2 <- exp(Wb)/(1+exp(Wb))      # 201*1 vector of P

Wb.sim <- W %*%t(b.sim)
P2.sim <- exp(Wb.sim)/(1+exp(Wb.sim))
dim(P2.sim)     #   61* 10000 simulated probabilities.
dim(apply(P2.sim , 1, sort))
dim(t(apply(P2.sim , 1, sort)))# 61*10000

P2.sim.graph <- t(apply(P2.sim,1,sort)) # 201*1000 simulated probabilities, sorted to identify the condidence bounds.


#
# Plot 
#


setwd("C:/Users/Jun/Dropbox/Elite Purge Project")

windows()
par(mar=c(5,5,4,1))
# A numerical vector of the form c(bottom, left, top, right) which gives the number of lines of margin to be specified on the four sides of the plot.
# The default is c(5, 4, 4, 2) + 0.1.
par(mfrow=c(1,2))


plot(tenure.range,P,type="l",lwd=3,col="black",  xlim=c(0, 11), ylim=c(0, 0.027) ,main="Probability of Coup Replacement", xaxt='n',ylab="Probability of Coup Replacement", xlab="Tenure")
lines(tenure.range,P2,type="l",lwd=3,col="black", lty=3)
axis(1, at=c(0,5,10,15), label=c(1, 6, 11, 16), col.axis="black")

legend(5, 0.027, c("Coup-Entry Leader", "No Coup-Entry Leader"),lwd=c(3,3), lty=c(1,3) , col="black")


# First-Difference Graphs

#windows()
#par(mfrow=c(1,1))

Diff.P.sim <- P.sim - P2.sim
dim(Diff.P.sim)  # 61*1000
rowMeans(Diff.P.sim)
dim(apply(Diff.P.sim,1,sort)) # 1000*61
dim(t(apply(Diff.P.sim,1,sort)))# 61*1000

plot(tenure.range, rowMeans(Diff.P.sim) ,type="l", xlim=c(0, 11),main="Effect of Coup Entry", ylim=c(-0.04, 0.05),xaxt='n', lwd=3,ylab="Pr(Coup Replacement|Coup-Entry)-Pr(Coup Replacement|No Coup-Entry)", xlab="Tenure")
abline(h=0)
lines(tenure.range, t(apply(Diff.P.sim,1,sort))[,9750],type="l",lty=3,lwd=2)
lines(tenure.range,  t(apply(Diff.P.sim,1,sort))[,250],type="l",lty=3,lwd=2)
#legend(x=500, y=0.6, bty="n", legend="Incumbency Advantage", lty=1, lwd=3, seg.len=0.75, col="blue", x.intersp=0.5)
axis(1, at=c(0,5,10,15), label=c(1, 6, 11, 16), col.axis="black")

# mtext('Simulation (a.k.a. CLARIFY)', outer=T, line=-2)
legend(5, 0.049, c("Effect of Coup-Entry", "95% CI"),lwd=c(3,2), lty=c(1,3) , col="black")

# save

savePlot(filename="Paper/purge_01012015_H2", type = "pdf")

 



##
## Figure 2 
##
 

b.sim <- read.table("C:/Users/Jun/Dropbox/Elite Purge Project/Data/VCVmatrix_01012015.csv", header=TRUE, sep=",")  # this is a 10000 draws of simulated coefficients.
dim(b.sim) # 10000 * 10


f <- purge4 ~  logmilex + IS_WAR  + lrgdpch + mildic + royal + party +  failedcoup2   + purge4yrs + I(purge4yrs ^2) + I(purge4yrs ^3) + coupentry2*log(tenure+1)
model1 <- glm(f, family=binomial(link = "logit"), data=dat)

coefficients(model1)
tenure.range <- seq(0, 12, 0.25)

# Coup-entry leader

Z <- cbind(1,   mean(dat$logmilex, na.rm=TRUE) ,  median(dat$IS_WAR , na.rm=TRUE), mean(dat$lrgdpch, na.rm=TRUE) ,   median(dat$mildic, na.rm=TRUE),
median(dat$royal, na.rm=TRUE), median(dat$party, na.rm=TRUE), median(dat$failedcoup2, na.rm=TRUE), mean(dat$purge4yrs, na.rm=TRUE),
 (mean(dat$purge4yrs, na.rm=TRUE))^2,  (mean(dat$purge4yrs, na.rm=TRUE))^3,  1,  log(tenure.range +1), 1*log(tenure.range +1))
coefficients(model1)
Zb <- Z%*%coefficients(model1)
P <- exp(Zb)/(1+exp(Zb))      # 201*1 vector of P

Zb.sim <- Z %*% t(b.sim)
P.sim <- exp(Zb.sim)/(1+exp(Zb.sim))   # 201*1000 simulated probabilities.
P.sim.graph <- t(apply(P.sim,1,sort)) # 201*1000 simulated probabilities, sorted to identify the confidence bounds.


# No coup-entry leader

W <- cbind(1,  mean(dat$logmilex, na.rm=TRUE) , median(dat$IS_WAR , na.rm=TRUE), mean(dat$lrgdpch, na.rm=TRUE) ,  median(dat$mildic, na.rm=TRUE),
median(dat$royal, na.rm=TRUE), median(dat$party, na.rm=TRUE), median(dat$failedcoup2, na.rm=TRUE), mean(dat$purge4yrs, na.rm=TRUE),   (mean(dat$purge4yrs, na.rm=TRUE))^2,  (mean(dat$purge4yrs, na.rm=TRUE))^3,
 0,  log(tenure.range +1), 0*log(tenure.range +1))

Wb <- W %*%coefficients(model1)

P2 <- exp(Wb)/(1+exp(Wb))      # 201*1 vector of P

Wb.sim <- W %*%t(b.sim)
P2.sim <- exp(Wb.sim)/(1+exp(Wb.sim))
dim(P2.sim)     #   61* 10000 simulated probabilities.
dim(apply(P2.sim , 1, sort))
dim(t(apply(P2.sim , 1, sort)))# 61*10000

P2.sim.graph <- t(apply(P2.sim,1,sort)) # 201*1000 simulated probabilities, sorted to identify the condidence bounds.





#
# Graph
#
windows()
par(mar=c(5,5,4,3))
# A numerical vector of the form c(bottom, left, top, right) which gives the number of lines of margin to be specified on the four sides of the plot.
# The default is c(5, 4, 4, 2) + 0.1.
par(mfrow=c(1,2))


plot(tenure.range,P,type="l",lwd=3,col="black", ylim=c(0.02, 0.13), xlim=c(0, 11), main="Probability of Military Purge", xaxt='n',ylab="Probability of Military Purge", xlab="Tenure")
lines(tenure.range,P2,type="l",lwd=3, lty=3, col="black")
axis(1, at=c(0,5,10,15), label=c(1, 6, 11, 16), col.axis="black")

legend(3, 0.125, c("Coup-Entry Leader", "No Coup-Entry Leader"),lwd=3, lty=c(1,3) , col="black")


# First-Difference Graphs

#windows()
#par(mfrow=c(1,1))

Diff.P.sim <- P.sim - P2.sim
dim(Diff.P.sim)  # 61*1000
rowMeans(Diff.P.sim)
dim(apply(Diff.P.sim,1,sort)) # 1000*61
dim(t(apply(Diff.P.sim,1,sort)))# 61*1000

plot(tenure.range, rowMeans(Diff.P.sim) ,type="l", ylim=c(-0.05, 0.2),xlim=c(0, 11), lwd=3,xaxt='n', main="Effect of Coup Entry", ylab="Pr(Purge|Coup-Entry)- Pr(Purge|Non Coup-Entry)", xlab="Tenure")
abline(h=0)
lines(tenure.range, t(apply(Diff.P.sim,1,sort))[,9750],type="l",lty=3,lwd=2)
lines(tenure.range,  t(apply(Diff.P.sim,1,sort))[,250],type="l",lty=3,lwd=2)
axis(1, at=c(0,5,10,15), label=c(1, 6, 11, 16), col.axis="black")

legend(4, 0.18, c("Effect of Coup-Entry", "95 % CI"),lwd=c(3,2), lty=c(1,3) , col="black")
#mtext('Simulation (a.k.a. CLARIFY)', outer=T, line=-2)
 

# save
savePlot(filename="Paper/purge_01012015_H3", type = "pdf")








##
## Appendix F 
##
 
#
# Figure 1  

b.sim <- read.table("Data/VCVmatrix_1_polity.csv", header=TRUE, sep=",")  # this is a 10000 draws of simulated coefficients.
dim(b.sim) # 10000 * 10

f <- coupexit~   lrgdpch + logmilex  + mildic + royal + party + failedcoup2_1  + coupentry2*log(tenure+1)

dat1 <- dat[dat$polity < 6 & dat$polity!=-88 & dat$polity!=-77 & dat$polity!=-66,]

model1 <- glm(f, family=binomial(link = "logit"), data=dat1)


#b.sim <- mvrnorm(10000,coefficients(model1),vcov(model1)) # 1000*19 matrix of simulated cofficients.
#dim(b.sim) # [1] 1000    9

coefficients(model1)
tenure.range <- seq(0, 12, 0.25)

# Coup-entry leader

Z <- cbind(1,   mean(dat$lrgdpch, na.rm=TRUE), mean(dat$logmilex, na.rm=TRUE),  median(dat$mildic, na.rm=TRUE),
median(dat$royal, na.rm=TRUE), median(dat$party, na.rm=TRUE),0,  1,  log(tenure.range +1), 1*log(tenure.range +1))
coefficients(model1)
Zb <- Z%*%coefficients(model1)
P <- exp(Zb)/(1+exp(Zb))      # 201*1 vector of P

Zb.sim <- Z %*% t(b.sim)
P.sim <- exp(Zb.sim)/(1+exp(Zb.sim))   # 201*1000 simulated probabilities.
P.sim.graph <- t(apply(P.sim,1,sort)) # 201*1000 simulated probabilities, sorted to identify the confidence bounds.


# No coup-entry leader

W <- cbind(1,      mean(dat$lrgdpch, na.rm=TRUE), mean(dat$logmilex, na.rm=TRUE),  median(dat$mildic, na.rm=TRUE),
median(dat$royal, na.rm=TRUE), median(dat$party, na.rm=TRUE),0,  0,  log(tenure.range +1), 0*log(tenure.range +1))

Wb <- W %*%coefficients(model1)

P2 <- exp(Wb)/(1+exp(Wb))      # 201*1 vector of P

Wb.sim <- W %*%t(b.sim)
P2.sim <- exp(Wb.sim)/(1+exp(Wb.sim))
dim(P2.sim)     #   61* 10000 simulated probabilities.
dim(apply(P2.sim , 1, sort))
dim(t(apply(P2.sim , 1, sort)))# 61*10000

P2.sim.graph <- t(apply(P2.sim,1,sort)) # 201*1000 simulated probabilities, sorted to identify the condidence bounds.


#
# Graph   
#


setwd("C:/Users/Jun/Dropbox/Elite Purge Project")

windows()
par(mar=c(5,5,4,1))
# A numerical vector of the form c(bottom, left, top, right) which gives the number of lines of margin to be specified on the four sides of the plot.
# The default is c(5, 4, 4, 2) + 0.1.
par(mfrow=c(1,2))


plot(tenure.range,P,type="l",lwd=3,col="black",  xlim=c(0, 9), ylim=c(0.007, 0.038) ,main="Probability of Coup Replacement (Polity)", xaxt='n',ylab="Probability of Coup Replacement", xlab="Tenure")
lines(tenure.range,P2,type="l",lwd=3,col="black", lty=3)
axis(1, at=c(0,4,8,10), label=c(1, 5, 9, 11), col.axis="black")

legend(4 , 0.038, c("Coup-Entry Leader", "No Coup-Entry Leader"),lwd=c(3,3), lty=c(1,3) , col="black")



# First-Difference Graphs

#windows()
#par(mfrow=c(1,1))

Diff.P.sim <- P.sim - P2.sim
dim(Diff.P.sim)  # 61*1000
rowMeans(Diff.P.sim)
dim(apply(Diff.P.sim,1,sort)) # 1000*61
dim(t(apply(Diff.P.sim,1,sort)))# 61*1000

plot(tenure.range, rowMeans(Diff.P.sim) ,type="l", xlim=c(0, 9),main="Effect of Coup Entry (Polity)", ylim=c(-0.06, 0.045),xaxt='n', lwd=3,ylab="Pr(Coup Replacement|Coup-Entry)-Pr(Coup Replacement|No Coup-Entry)", xlab="Tenure")
abline(h=0)
lines(tenure.range, t(apply(Diff.P.sim,1,sort))[,9750],type="l",lty=3,lwd=2)
lines(tenure.range,  t(apply(Diff.P.sim,1,sort))[,250],type="l",lty=3,lwd=2)
#legend(x=500, y=0.6, bty="n", legend="Incumbency Advantage", lty=1, lwd=3, seg.len=0.75, col="blue", x.intersp=0.5)
axis(1, at=c(0,4,8,10), label=c(1, 5, 9, 11), col.axis="black")

# mtext('Simulation (a.k.a. CLARIFY)', outer=T, line=-2)
legend(4, -0.03, c("Effect of Coup-Entry", "95% CI"),lwd=c(3,2), lty=c(1,3) , col="black")


# save

savePlot(filename="Paper/purge_01012015_H2_polity", type = "pdf")




##
## Figure 2 
##
 

b.sim <- read.table("C:/Users/Jun/Dropbox/Elite Purge Project/Data/VCVmatrix_01012015_polity.csv", header=TRUE, sep=",")  # this is a 10000 draws of simulated coefficients.
dim(b.sim) # 10000 * 10


f <- purge4 ~  logmilex + IS_WAR  + lrgdpch + mildic + royal + party +  failedcoup2   + purge4yrs + I(purge4yrs ^2) + I(purge4yrs ^3) + coupentry2*log(tenure+1)

dat1 <- dat[dat$polity < 6 & dat$polity!=-88 & dat$polity!=-77 & dat$polity!=-66,]

model1 <- glm(f, family=binomial(link = "logit"), data=dat1)

coefficients(model1)
tenure.range <- seq(0, 12, 0.25)

# Coup-entry leader

Z <- cbind(1,   mean(dat$logmilex, na.rm=TRUE) ,  median(dat$IS_WAR , na.rm=TRUE), mean(dat$lrgdpch, na.rm=TRUE) ,   median(dat$mildic, na.rm=TRUE),
median(dat$royal, na.rm=TRUE), median(dat$party, na.rm=TRUE), median(dat$failedcoup2, na.rm=TRUE), mean(dat$purge4yrs, na.rm=TRUE),
 (mean(dat$purge4yrs, na.rm=TRUE))^2,  (mean(dat$purge4yrs, na.rm=TRUE))^3,  1,  log(tenure.range +1), 1*log(tenure.range +1))
coefficients(model1)
Zb <- Z%*%coefficients(model1)
P <- exp(Zb)/(1+exp(Zb))      # 201*1 vector of P

Zb.sim <- Z %*% t(b.sim)
P.sim <- exp(Zb.sim)/(1+exp(Zb.sim))   # 201*1000 simulated probabilities.
P.sim.graph <- t(apply(P.sim,1,sort)) # 201*1000 simulated probabilities, sorted to identify the confidence bounds.


# No coup-entry leader

W <- cbind(1,  mean(dat$logmilex, na.rm=TRUE) , median(dat$IS_WAR , na.rm=TRUE), mean(dat$lrgdpch, na.rm=TRUE) ,  median(dat$mildic, na.rm=TRUE),
median(dat$royal, na.rm=TRUE), median(dat$party, na.rm=TRUE), median(dat$failedcoup2, na.rm=TRUE), mean(dat$purge4yrs, na.rm=TRUE),   (mean(dat$purge4yrs, na.rm=TRUE))^2,  (mean(dat$purge4yrs, na.rm=TRUE))^3,
 0,  log(tenure.range +1), 0*log(tenure.range +1))

Wb <- W %*%coefficients(model1)

P2 <- exp(Wb)/(1+exp(Wb))      # 201*1 vector of P

Wb.sim <- W %*%t(b.sim)
P2.sim <- exp(Wb.sim)/(1+exp(Wb.sim))
dim(P2.sim)     #   61* 10000 simulated probabilities.
dim(apply(P2.sim , 1, sort))
dim(t(apply(P2.sim , 1, sort)))# 61*10000

P2.sim.graph <- t(apply(P2.sim,1,sort)) # 201*1000 simulated probabilities, sorted to identify the condidence bounds.





#
# Graph
#
windows()
par(mar=c(5,5,4,3))
# A numerical vector of the form c(bottom, left, top, right) which gives the number of lines of margin to be specified on the four sides of the plot.
# The default is c(5, 4, 4, 2) + 0.1.
par(mfrow=c(1,2))


plot(tenure.range,P,type="l",lwd=3,col="black", ylim=c(0.02, 0.13), xlim=c(0, 11), main="Probability of Military Purge (Polity)", xaxt='n',ylab="Probability of Military Purge", xlab="Tenure")
lines(tenure.range,P2,type="l",lwd=3, lty=3, col="black")
axis(1, at=c(0,5,10,15), label=c(1, 6, 11, 16), col.axis="black")

legend(3, 0.125, c("Coup-Entry Leader", "No Coup-Entry Leader"),lwd=3, lty=c(1,3) , col="black")

# save
#savePlot(filename="Dissertation_slide_H3", type = "pdf")


# First-Difference Graphs

#windows()
#par(mfrow=c(1,1))

Diff.P.sim <- P.sim - P2.sim
dim(Diff.P.sim)  # 61*1000
rowMeans(Diff.P.sim)
dim(apply(Diff.P.sim,1,sort)) # 1000*61
dim(t(apply(Diff.P.sim,1,sort)))# 61*1000

plot(tenure.range, rowMeans(Diff.P.sim) ,type="l", ylim=c(-0.05, 0.2),xlim=c(0, 11), lwd=3,xaxt='n', main="Effect of Coup Entry (Polity)", ylab="Pr(Purge|Coup-Entry)- Pr(Purge|Non Coup-Entry)", xlab="Tenure")
abline(h=0)
lines(tenure.range, t(apply(Diff.P.sim,1,sort))[,9750],type="l",lty=3,lwd=2)
lines(tenure.range,  t(apply(Diff.P.sim,1,sort))[,250],type="l",lty=3,lwd=2)
axis(1, at=c(0,5,10,15), label=c(1, 6, 11, 16), col.axis="black")

legend(4, 0.18, c("Effect of Coup-Entry", "95 % CI"),lwd=c(3,2), lty=c(1,3) , col="black")
#mtext('Simulation (a.k.a. CLARIFY)', outer=T, line=-2)

# save
#savePlot(filename="Dissertation_12172012_H2_2_difference", type = "pdf")


# save
savePlot(filename="Paper/purge_01012015_H3_polity", type = "pdf")







##
## End ; 
##


